arXiv: 1506.08517vl [math.NA] 29 Jun 2015 


A GENERALIZATION OE THE DIVIDE AND CONQUER ALGORITHM FOR THE 
SYMMETRIC TRIDIAGONAL EIGENPROBLEM 

DO YOUNG KWAK* AND JAEYEON KIM t 


Abstract. In this paper, we present a generalized Cuppen’s divide-and-conquer algorithm for the symmetric 
tridiagonal eigenproblem. We extend the Cuppen’s work to the rank two modifications of the form A = T + /3iwiw^ + 
^ 2 W 2 wJ, where T is a block tridiagonal matrix having three blocks. We introduce a new deflation technique and 
obtain a secular equation, for which the distribution of eigenvalues is nontrivial. We present a way to count the number 
of eigenvalues in each subinterval. It turns out that each subinterval contains either none, one or two eigenvalues. 
Furthermore, computing eigenvectors preserving the orthogonality are also suggested. Some numerical results, showing 
our algorithm can calculate the eigenvalue twice as fast as the Cuppen’s divide-and-conquer algorithm, are included. 


Key words. Eigenvalues of symmetric tridiagonal matrix, Cuppen’s divide and conquer method, 
Rank two modifications. Secular equation. 
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1. Introduction. Given a n x n tridiagonal real symmetric matrix A, the symmetric tridiagonal 
eigenproblem is to find all the eigenvalues and the corresponding eigenvectors of A. In this paper 
we generalize the Cuppen’s divide and conquer algorithm (CDC) [1] for the symmetric tridiagonal 
eigenproblem by considering a rank two modification. The CDC algorithm starts with decomposing 
A as a rank one modification 


A = 


r Ti 

Q 

OI 

_i 

T2 


+ /3ww^, 


w — Ofc + efc_|_i, 


( 1 . 1 ) 


where Ti,T 2 are symmetric tridiagonal matrices of order k^l, n — k^l respectively, and /3 7 ^ 0 
is (k, fc -b 1) entries of A. By diagonalizing Ti as Ti = QiDiQj and deflating the cases when the 
eigenvalues of A coincide with the diagonals of Di, secular equation of CDC is derived. By using the 
eigenvalues obtained from the secular equation, we can compute eigenvectors. 

It is tempting to generalize the CDC algorithm to the following rank two modifications: 

A = T +/3iwiwf-b/32W2 wJ, (1.2) 

where T is block tridiagonal, and wi = -b W 2 = Then we mimic the CDC 

algorithm to compute eigenvalues and eigenvectors of (O by using spectral decomposition of T^’s. 
We deflate the cases when some of the eigenvalues of A coincide with the diagonals of Di. Then we 
can derive the secular equation. By computing the eigenvalues from the secular equation, we can 
finally obtain the eigenvectors of A. 

2. Dividing Step. We rewrite A blockwise as follows: 



OI 

OI 


ill I 0 


1 

o 

lO 

lO 

A = 

OI 

OI 

+ /3i 

^ F I Q 

+ /32 

Q : ill 


lO 

lO 


lO 

lO 

lO 


Q I ^ I 1 


where Ti, T 2 and T 3 are symmetric tridiagonal matrices of order fci ^ 1 , ^2 ^ 1 , fca ^ 1 respectively, 
with -b ^2 + fca = n and 7 ^ 0, /32 7 ^ 0 are (fci, fci -b 1), (^ 2,^2 + 1) entries of A respectively. The 
computation of the spectral decomposition of each Ti 


Ti = QiDiQj, i = l,2,3 
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is cheaper than (about 2/3) those of CDC where Di are diagonal matrices and QjQi = /. Let 


■ Qi I Q I 0 ■ 


■ Di I 0 I 0 ■ 

■■QTQ 2 TQ" 

II 

■■■q'T'dVT'q’" 

■■qTttq;" 




Since Q is an orthogonal matrix, eigenvalues of AQ are same as those of A. We have 

AQ = Q^{T + /3iwiwf + /32W2w|’)Q = D + /3iVivf + 132^2^2, (2.2) 


where Vi = Q^Wi, V 2 = Q^W 2 . Then 


( last column of Qj\ / Q \ 

first column of Q 2 , V 2 = last column of Q 2 

0 ) yfirst column of Qg ) 

Let X be an eigenvector of Q^AQ corresponding to an eigenvalue A. Then we have 

{D + /3ivivf + /32V2 v|’)x = Ax, 

which is equivalent to 


(2.3) 


{D - A/)x =-Pi (vf x)Vi - p 2 (vf x)V 2 . 


(2.4) 


If {D — XI) is invertible, x can be expressed as a linear combination of Ui and U 2 where Ui = 
{XI — D)~^'Vi, U 2 = {XI — D)~^V 2 . Using this expression of x, we can derive the secular equation 
(see Section S]). 

Note that the matrix {D — XI) is singular if and only if the eigenvalues of A coincide with the 
diagonals oi D. In order to exclude those cases, we deflate them by computing the corresponding 
eigenvectors. This will be explained in detail in the next section. 

3. Deflation. In this section we will determine the cases when some of the diagonal entries of 
D are eigenvalues of A. Expressing equation (12.41) componentwise, we see 


X)xi^ 


//3i(vfx)?;ii\ 


/ P2{'V2^)v2i\ 

X)X2 

+ 

Pl{vjx)vi2 

+ 

P2{v2^)v22 

X)XnJ 


\Pl{vJ :x.)vin J 


\P2{'V2 ^)v2n J 


= 0 . 


(3.1) 


We will check row by row whether each diagonal entry is an eigenvalue or not. We first assume 
di is an eigenvalue. Based on the values of wii,U 2 i, we can divide it into four cases as follows: 

1. Assume vu = 0 and V 2 i = 0. Then is an eigenvector for di. 

2. Assume vu = 0 and V 2 i 0. Since P 2 ^ 0, we have v|’x = 0. Hence, deleting the Uth row in 
equation (EB, we obtain the following equation. 

{D' - A/)x' + /3i(vl(vl)^)x' = 0, (3.2) 

where iA',Vg,x' are obtained by omitting the i-th entry. In this case we can proceed exactly 
as in Cuppen’s divide and conquer method to check whether di is an eigenvalue or not. This 
algorithm is summarized as follows. [I] 

(a) If there exists j such that dj = di,vij = 0, then di is an eigenvalue of D' (and hence of 
D) and x' = e' where e' is obtained from by omitting the i-th entry. 
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(b) If there exists j, k such that dj = dk = di, then di is an eigenvalue and 


/ 

■ 0 ■ 

\ 


-Vlk 

jth 


Vij 

kth 

V 

0 

/ 


where / is obtained by omitting the i-th entry. 

(c) If there is no j such that dj = di, we obtain the secular equation of (13.2|) from the 
Cuppen’s divide and conquer method as follows: 


/ 11 (A)1 - /3i 

q^i 


V 


2 

I 9 


{x-d,y 


(3.3) 


If fii{di) = 0, then di is an eigenvalue and a;(. = for all r ^ i. Hence we define 

/ 11 (A) as a discriminant for di. 

For all three cases, we can find Xi by solving v^x = 0 with given x' above. Therefore, di is 
an eigenvalue and the corresponding eigenvector is x. 

3. Assume vu ^ 0 and V 2 i = 0. We can apply the same method as above and obtain the following 
a discriminant: 


/ 12 (A) -I-P 2 Y. 

q^i 


V 


2 

2 g 


(A-d,)- 


(3.4) 


4. Finally assume vu ^ 0 and V 2 i ^ 0. Suppose the diagonal entry di is repeated exactly p 
times. Rearranging the indices, let us assume dmi = dm 2 = ■ ■ ■ = dmp- From equation (13.IL 
we see the j-th row, for j = mi, ■ ■ ■ , mp, satisfy the following equations: 


{ /3l(vfx)uimi +/32(v|’x)u2mi = 0 

/3l(vfx)uimp +/32(vfx)u2mp = 0. 


(3.5) 


In matrix form, this can be written as 


Cz = 0, 


where 


C := 


Vlmi 

yim2 


V2mi 

'1^27712 


is p X 2 and z 


Vlrrip '^2mp 


/3i(vfx) 


is 2 X 1 vector. 


Depending on the rank of C, these are divided into two cases: 
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(a) Assume rank{C) = 2. Then equation (13.5|) has the trivial solution, which means v^x = 0 
and vi’x = 0. Since {dq — di) ^ 0, we see from equation (13.11) that Xq = 0, for all 
q 7 ^ nil , TO 2 ... nrip. Therefore, 

/ X — ^Im2^m2 Vlmp^nip — 0 13 61 

[^^2 X — X2mi^mi T ‘^2m2^m2 V2mp^mp — 0- 

From these equations, we can find exactly p — 2 independent eigenvectors. In this case 
we still have two identical diagonal entries which are not deflated. 

(b) Assume rank{C) = 1. Consider the case when vfx = 0. In this case it follows that 
^2 X = 0 and we can find exactly p — 1 nontrivial solutions associated with eigenvalue 
di- Since Xg = 0 for every q ^ mi, m 2 ,..., mp, by the same way as above we can deflate 
the corresponding p — 1 rows and columns. After the deflation, we would have only one 
diagonal entry of D which is equal to di . Without loss of generality, we may assume that 
this diagonal entry is dmi- 

Next consider the case when v^x 7 ^ 0. In this case, it follows that vi’x 7 ^ 0. Fix any 
such X. Then we have 


v^x = I, 


(3.7) 


for some 1^0. From equation dSSl), we obtain 

/32W2mi(vi’x) / 32 IV 2 n 


Vi X =- 


Substitute this value of x into equation (I33D, we get 

''9 ~ 


( )^l9 + ^2^t’2g ^ VirmV2q - V2rm 


Vlq 


dm\ dq 


^Imi 


dmi dq 


for all q 7 ^ mi. Solving equation (j3.8l) using these Xq, we obtain Xm^ as 


(3.8) 


(3.9) 


q^mi 

_ I32lv2m.i ^2l '^2m\'^lq 


q^mi 


^mi dq 


(3.10) 


Substituting the values Xq in p.9l) and p.IOl) for all g = I, • ■ • , into the equation (13.7p . 


I =v^x 


— ^ ^ '^2q^q ~\~ '^2mi^mi 

q^mi 


hi 


- E ^29 




q^mi 


dm I dq 


l32lV2mi 


!« 2 mi 

?:::r ^ 


Vlq 

q^mi 
,,2 ^,2 


'^lmi'^2q '^2mi'^lq 

dmi dq 


i 82 *t' 2 mi hi Ulmil^iq 2,VimiV2miVlqV2q + V2m,^vfq 


q^mi 


dmi ^q 
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Since Z ^ 0, we obtain 

P2V2mi 1^2 


1 




E 


- ‘2-Vlm^V2rmVlqV2q + 
dnii dq 


= 0 . 


(j/mi 

Therefore, we define the following function as the discriminant in this case: 

MX) := + &»L. - ftA i: 


q^mi 


If / 2 (di) = 0, then di is an eigenvalue and the corresponding eigenvector is x given in 
(1531) and (13.101) . Otherwise, di cannot be an eigenvalue. 

Since we have gone through every possible cases, di cannot be an eigenvalue after deflation. 

4. Secular Equation. We now assume that we have deflated all the cases as above. 

Theorem 4.1. The eigenvalues of 


D + /3ivivf + /32V2V^ 


are the roots of the secular equation defined by 


/(A):= 


q—1 r—q-\-l 


{viqV2r “ VirV2q)'^ 

{X-dg){X-dr) 


E 

9=1 




1? 


■ I32v'^ 


2q 


X dn 


1 . 


(4.1) 


Proof Since the di is not an eigenvalue for all i, {D — XI) is invertible. From (12.4L we have 

x= (A/- i:>)“^(/3i(vfx)vi +/32 (v^x)v 2) (4.2) 

:=aui+6u2, (4.3) 

where a = /li(vfx), b = /32(vjx), Ui = {XI — U2 = {XI — D)~^V 2 . Substituting (14.31) into 

(131) . we obtain 

{—a + a/3ici + 6/3 iC 3 )vi + {—b + 6/32C2 + a/32C3)v2 = 0, (4.4) 

where ci = vfui, C2 = vju2, C3 = vfu2 = vjui. 

If Vi is multiple of V2, this problem is reduced to the original divide and conquer method of 
Cuppen. Therefore, we assume that Vi is not multiple of V2. Hence we have 

— a + a/3ici + &/3 iC 3 =0, —b + b/32C2 + a(32C3 = 0. (4.5) 

Since a(l — /3iCi) = b/Sic^ and 6(1 — /32C2) = aj32C3, by multiplying we obtain 

a6(l - /3ici)(l - /32C2) = a6/3iC3/32C3. (4.6) 

First we assume ab 7^ 0. Then we can cancel ab and we obtain 

(-1 +/3ici)(-l +/32C2) =/31/32C3. (4.7) 

Here, we define the secular equation /(A) as follows: 

/(A) := /3i/32(ciC2 - C 3) - /3iCi - /32C2 + 1. 


(4.8) 
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By substituting values of ci, C 2 and C 3 , we obtain 


n 2 ^ 2 * 


9=1 


9=1 


^9=1 


9=1 


A dn 


(4.9) 


=AAi: 1 : 


{viqV2r “ 'l’lrl’2g)^ /^l^^lg + P^vlq 

^ A-d„ 


9=lr=9+l ,= 1 “9 


If a = 0, & 7 ^ 0, we see C 3 = 0 and 1 — / 32 C 2 =0. If a 7 ^ 0 ,6 = 0, we see C 3 = 0 and I — /3ici = 0. We 
see these cases also satisfy the secular equation (HU). □ We note that the derivative of /(A) is 

{2\ - dq - dr){viqV2r - VlrV2qY l^lVlq + l32V2q 

(A - dqY{X - drf ^ i^-dq^ ■ 


/'(A) = -/3i/32^ ^ 

q—1 r—q-\-l 


5 . Computing eigenvalues from secular equation. Since triple or more identical diagonal 
entries are deflated by Section 3, cases (1), (2)-(b), (3)-(b), (4)-(a),(b), there can be at most two 
identical diagonal entries for each diagonal value of D. Fix an index i = iq. If there exists an index 
jo ^ *0 such that djq = di^, we call this diq multiple. Otherwise, we call dig simple. 

First, we will check the sign of lim /(A) and lim /(A) which are determined by the sign of/(di^+e) 

*0 *0 _ 

for small e. By checking the dominating term in (14.8L we see it has the same sign as 


P 1 P 2 Y + 1^2.^210 

^ ^ -dq + e)e e 

9#*o 


dip simple 


dig multiple. 


For multiple dig, we know that {vijgV 2 ig — viigV 2 jo)'^ 7 ^ 0 from (2)-(a), (3)-(a), (4)-(b) of Section 3. 
Hence the sign of lim /(A) and lim /(A) will be determined by the coefficients of - and which 

X^df A->-dr 

*0 ^0 


are given as follows: 


at ■= 


-P 1 P 2 Y. 


P 1 P 2 , 


9=1 

95^*0 


(^l9^^2yo - Vl,ioV2q)‘^ 

dq — dig 


- PlVUg - P 2 V: 


2 

2io ’ 


dig simple 


dig multiple. 


(5.1) 


9ig 


P 1 P 2 


P 1 P 2 , 


^ {viqV2,io - Vl,ioV2q f 

9=1 
9#*0 


dq dig 


+ /3lOio + ^ 211 ; 


2ig ^ 


dig simple 
dig multiple. 


(5.2) 


Now, we are going to determine the number of roots for each interval {di,di+i). First, we are 
going to deal with the cases without multiple diagonal entries ('Section [54^ . After that, we are going 
to deal with the cases having some multiple diagonal entries ('Section lOl) . 
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5.1. Cases without multiple diagonal entries. Suppose that (n — m) diagonal entries are 
deflated in the steps (1) - (4) of Section 3. For the simplicity of presentation, let us sort and reindex 
the remaining m diagonal entries in an ascending order as follows: 


0 = do < di < ■ ■ ■ < dm < dm+i = C)0. 

We let Iq := {dq, dq+i) and define Iqr as the union of (r — q) intervals as 

r—1 

Iqr • — li — (dg , dg-|_i) U (dg-|_i, dq-\-2 ) U ■ * ■ U (^dr — l , dr ), 

i=q 

where q = 0,1,2, ,m and r = q + 1,... ,'m + l. In addition, we define the following matrices of sizes 
fci + k 2 and fc 2 + ko respectively. 


s,; := 


Dr 


D 


i+l 


+ /3iV-v' , i = l,2 


where v' is obtained by omitting zero vector part from (12.31) . Then we see 

Di I 


A = 


\ 



Do \ 


+ ^2V2V2 = 


Bo 


T^iViVi . 


(5.3) 


(5.4) 


We assume the following spectral decomposition of Bii 

Br=Q'D'Q^, Q'^Q'=I, 1 = 1,2 


and define 


\ Q[\ ] 


' I \ 

i I 

5 J Cj 

1 Q 2 


Since and i ?2 are orthogonal matrices, A has same eigenvalue as R"[ARi and RjAR 2 : 


R^ARi 



+ P2Z2Z2, 


R 2 AR 2 



+ /3izizf 


where Z 2 = i?fv 2 ,Zi = Also, we define the interval /' for j = 1, 2, • • • ,m— las 


I':= (d",d^,) 


(5.5) 


where d^ are diagonal entries from and Do (or Di and H^) in an ascending order. Since the 
matrices in (15.51) are in the rank one modification form, we see from the Cuppen’s divide and conquer 
algorithm that there exists exactly one eigenvalue of Rf ARi (and hence of A) in each . 



Fig. 5.1. /fc most two eigenvalues when d'J is in (left) and has at most one eigenvalue when d^ is 

not in Ik (right) 
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Lemma 5.1. Interval Ijj+k can have at least k — 1 roots and at most fc + 1 roots for j = 
1, 2, • • • , m — 1 and fc = 1, 2, • • ■ — j. 

Proof. Fix the index j = jo and fc = fco so that the Ijgjg+ko has fco + 1 diagonals including both 
ends. If djg+ko is from Di or D 2 , choose Bi = Bi. Otherwise, choose Bi = B 2 . Assume that Bi has 
I diagonal entries. Depending on whether djg is from those I diagonals or not, we can divide in into 
two cases. 

Assume that djg is one of the diagonals of Bi. Then from (15.31) we conclude by GDC algorithm 
that Bi has exactly I — 1 eigenvalues in Ijg^g+kg- Since each of ^ — 1 eigenvalues is one of d”, 
there will be at least (Z — 1) + (fco + 1 — 1 ) such d"s in Ijgjg+ko- {ko + 1 — I d"’s are from D 3 
or Di). Then we have 

^s-i ^ cijg < dg < ■ ■ ■ < dg_^i.g_i < djg+ko < ^s+ko (^-h) 

for some s. Since each /' has exactly one eigenvalue of A for q = s — - ,s + fco — 1, 

Ijgjg+ko can have at least fco — 1 and at most fco + 1 eigenvalues of A (see Figure [SH]) 
Assume that djg is not one of the diagonals of Bi. (i.e., djg G D 3 or djg G Di.) Then Bi can 
have I — 1 01 I eigenvalues in Ijgjg^ko by CDC algorithm. Since each eigenvalues of Bi is one 
of d", there will be (Z — 1) + (fco + 1 — Z) or / + (fco + 1 — 1 ) such d"s in Ijgjg+ko including djg. 
Then we have 

r d" = djg < d"+i < • • • < <+feo_i < djg+kg < d"+fe^ if fco such d"s 
1 d” = djg < d;;+i < • • • < d"+fe^ < djg+kg < d"+fco+i if fco + 1 such d"s 

for some s. Then Ijg jgj.kg can have fco — 1 or fco eigenvalues of A for the first case and can 
have fco or fco + 1 eigenvalues of A for the second case. 

In both cases Ijgjg+ko fco ~ 1 or fco or fco + 1 eigenvalues of A, which completes the proof of 
Lemma. □ From this Lemma, we can say that there can be at most two eigenvalues of A in each 
interval Ij,j = l,--- ,m — I. Let us classify every intervals into two groups S~^ and S~ as follows: 

f Ijk G S+, if > 0 

\ Ijk G S-, if < 0 ■ 

Since every interval can have at most two roots by Lemma l5. 1 1 and the secular equation is continuous 
on every interval Ij, we can conclude that every interval in S~ has one root by intermediate value 
theorem. Clearly, every interval Ij in S'®" has two roots or no root. By using mathematical induction, 
we can prove the following lemma. 

Lemma 5.2. For any interval Ijj+k in S~ has exactly fc roots where j = l,2,...,m — I and 
k < m — j. Also, for any interval Ijj+k in has fc — I or fc + 1 roots where j = 1,2,... ,m — 1 and 
k < m — j. 

Proof. We have already seen that each interval Ij satishes the lemma. Next, fix j = jo and 
suppose that the statements of the lemma holds when fc = fco. 

(1) Assume /jojo+feo G S~ and Ijg+kg G S~. Since g^g+^g and has opposite sign, 

Ijgjg+kg+i G S~. By induction, there will be fco roots in Ijgjg+kg and one root in Ijg+kg. 
Therefore, Ijgjg+kg+i will have fco + 1 roots and fc = fco + I satishes the lemma. 

(2) Assume Ijgjg+kg G S~ and Ijg+kg G S®". Then we have Ijgjg+kg+i G S®". By induction, there 
will be fco roots in Ijgjg+kg and two or no roots in Ijg+kg. Therefore, Ijgjg+kg+i will have fco 
or fco + 2 roots and fc = fco + 1 satishes the lemma. 

(3) Assume Ijgjg+kg G 5®" and Ijg+kg G S~. Then we have Ijgjg+kg+i G 5®". By induction, there 
will be fco — I or fco + 1 roots in Ijgjg+kg and one root in Ijg+kg. Therefore, Ijgjg+kg+i will 
have fco or fco + 2 roots and fc = fco + 1 satishes the lemma. 


( 1 ) 


( 2 ) 


A generalization of the Cuppen’s algorithm 


9 


(4) Assume Ijojo+ko S 5'“'' and Ijg+ko G S^. Then we have Ijgjg+ko+i G S~. By induction, there 
will be fco — 1 or fco + 1 roots in Ijgja+ko and two or no roots in Ijg+ko- Then Ijgjg+ko+i can 
have fco — 1 or fco + 1 or /co + 3 roots. However, this interval only can have fco + 1 roots by 
Lemma 15.11 and fc = fco + 1 satisfies the lemma. 

Hence the statements of the lemma holds when fc = fco + 1. □ Since the secular equation can have at 
most n roots, from Lemma IQ we see that Iq U Im can have at most two roots. 

Theorem 5.3. Number of roots for each interval is determined by the sign of g~ for i = 
1,2 ,... ,m. 

Proof. For each interval /j, j = 1, 2,... m — 1, we can determine the exact number of roots by the 
following: If Ij G S~, there will be only one root in Ij. If Ij G S'^, there can be at most two eigenvalues 
by Lemma l5.Il We can determine the exact number of roots by evaluating the sign of /(Aq) where Ao 
is a solution of /'(A) = 0. If the sign of /(Aq) is the same as those of /(d^) and /(d~^i), there will 
be no root. On the other hand, if the sign of /(Aq) is the opposite of f{d'^) and f{dj_^_-^), there will 
be two roots in Ij. Therefore, we know the number of roots in Iim since Iim = /i U • ■ • Im-i- Now, 
we determine the number of roots in /q and Im. Depending on the sign of /(d^) and f{dm), we can 
divide it into four cases as follows: 

(1) If gi > 0, 5“ > 0, we see Iim G S~^ and 5+ < 0. Since limA^oo /(A) = 1, there exists one root 
in Im- Since Iim can have m — 2 or to roots by Lemma 15.21 Iim bas to — 2 roots. Therefore, 
Iq has one root. 

(2) If g^ < 0, > 0, we see hm G S~ and 5+ < 0. Since limA->.oo /(A) = 1, there exists one 

root in Im- Since Iim has to — 1 roots by Lemma 15.21 Iq has no root. 

(3) If gi > 0, < 0, we see hm G S~ and > 0. Again by the fact limA^oo /(A) = 1 and 

gm > 0, Im can have two or zero roots. Since hm has to — 1 roots by Lemma [5^ Im cannot 
have two roots. Hence Im has no root and Jo must have one root. 

(4) If gi < 0, < 0, we see hm G S'^ and 5+ > 0. Then there can be to — 2 or to roots in 

hm by Lemma 15.21 Since Im belongs to S~^, there will be two or no roots in Im- Hence, by 
counting, we see /q also has two or no roots. Depending on the sign of Pi and P 2 , this is 
classified into the following cases: 

(a) Assume /3i > 0 and P 2 > 0. If dm is from Di or D 2 , choose Bi = Bi. Otherwise, choose 
Bi = B 2 - Applying the result of CDC to the form (15.3L we see dm < for some r. It 
is well known that when we apply the CDC algorithm for (ED then there is a root at 
the left end interval when P is negative and there is a root at right end interval when P 
is positive. Applying this result of CDC to the form (15.31) . we see there is at least one 
root in the interval (d^dm+i). This root belongs to Im- Since Im can have two or no 
roots, we see that Im has 2 roots and hence Iim has to — 2 roots. 

(b) Assume Pi < 0 and P 2 < 0- Arguing exactly the same way as in case (a), we conclude 
that the interval Iq has at least one root. This implies that Jo has 2 roots and Iim has 
TO — 2 roots. 

(c) Assume P 1 P 2 < 0. Similarly, we can see that neither Iq nor Im has two roots. This 
implies that hm has to roots and there is no root in Iq and Im- 

□ 


5.2. Cases with multiple diagonal entries. For the multiple diagonal case, the results of 
previous section can be applied in a similar manner. Suppose that di = di+i for some i. From (15.11) 
and (j5.2p . we see the sign of gi^ and g~ are the same. 

Again by applying the idea used in the proof of Lemma 15.11 and 15.21 we see that there should be 
one more root in one of h-i.i or Ii+i,i+ 2 - 
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Typical case of (3) Typical case of (4)-(b) 

Fig. 5.2. Typical examples of Theorem El 


6 . Eigenvector. Since deflation steps were treated in Section 3, it suffices to calculate eigenvec¬ 
tors using equation (14.31) . From (14.51) we have 

/-l + /3ici ^1C3 0 (6 1) 

P2C3 -l + P2C2j{b) - 

Since this system has a nontrivial solution, the determinant of the matrix must be zero, which is 
nothing but (14.71) . Eigenvectors corresponding to simple eigenvalues can be obtained by substituting 
the eigenvalues into (14.51) to solve for the ratio of a to b. For the case of multiple(double) eigenvalues, 
two linearly independent solutions exist. Hence all the entries of the matrix vanish. So any two 
nonzero numbers a, b are solutions. Hence we can take Ui and U2 as the corresponding eigenvectors. 

In practice, the eigenvectors computed this way lose orthogonality as in the original Cuppen’s 
divide and conquer algorithm when eigenvalues of A are close to those of T^. To fix these problems, 
we try two methods: 

6.1. Method 1 - Calculating Vi,V 2 corresponding to the computed eigenvalues. This 
is a natural extension of Gu and Eisenstat [5] fixing the orthogonality problem to the case of rank two 
modifications. However, as it turns out, we run into the lack of equations to find Vi and V2. To see 
why, first we need the following Lemma which can be found in the standard textbook. 






























A generalization of the Cuppen’s algorithm 
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Lemma 6.1. We have 
(1) det{A += det{A){l-\-A 

T\-l _ _ §_A^^vv'^A~^ 


T A-l. 


( 2 ) {A + Pw 


l+/3v^A-iv • 

Suppose that we could find two vectors Vi, V2 such that {Ai} are the exact eigenvalues of the new 
rank-two modification matrix D + + /32V2v|’. Proceeding as in Section 3 of [5] we are led to 

investigate the term nr=i(^* ~ d)- Using the Lemma, we have 


]^(Ai - ^J)= det{D - fil + /3ivivf -k P 2 ^ 2 ^ 2 ) 

= det{D — fj.1 + ,5iVivf )(1 -I- /32v|’ {D - fil + /3ivivf )“^V 2 ) 

= det{D - /r/)(l -k ftvf {D - /rJ)“ Vi)(l -k ^ 2^2 {D - fj.1 + /3iVivf )“iv2) 
= det{D - + pivj{D - /rJ)-ivi) (l -t ^2^1 {(£> - ^/)-i - 


-}v.) 


= det{D — fj.1) 1^1 -I- 

n / n 

= - Ai) ( 1 +X! 


<?=i 

n n 


Pivjq + P 2 vlq 

dq- H 

Pivfq + P 2 V 2 q 

dq- 


n 2 n 2 / n \ 


1 




q—1 r—q-\-l 


{viqV2T - VlrV 2 qY 

{dq - n){dr - ^i) 


^ ^9 “ ^ / 


= W{di - IJ) +^W{d, - ^x){Pivl^ + P2vlq) + ^1^2^^^ X! n 


2 = 1 


q—1 2=1 
i^q 


q=l r=q+l i=l 
i^q,r 


Setting /r = dk, k = 1, • • • , n, we obtain relatively simple equations. However, unlike the rank one 
modification, it is obvious that we cannot solve these equations for Vi and V2- We tried some ad-hoc 
method to compute approximations of vi and V2. The advantage of this method is that it requires 
only O(n^) operations. Since this approximation leads the orthogonality problem, the results were 
not so satisfactory for n > 80. 

6.2. Method 2 - Repeated applications of Gu and Eisenstat. We apply the methods used 
in Gu and Eisenstat [5] by regarding our decomposition as a repeated rank one modifications. 

(1) Apply the Gu’s method to Bi of (15.31) (resp. B 2 ) then 

(2) apply the Gu’s method to the first expression of A in (15.4|) (resp. second expression of A). 
Although this method needs O(n^) operations, the orthogonality problem does not arise. Hence, we 
can apply this algorithm in a recursive way. 


7. Numerical results. In this section we compare the execution time and accuracy measures 
between GDC and our rank two modification divide and conquer(RTDC). The algorithm was run on 
a 32 bit laptop computer which has machine epsilon e = 2.2204 x 10“^®. As a test matrix, we have 
chosen a typical matrix arising from solving Laplace equations on squares using the finite difference 
method. The test matrix is 



' B 

1 

10 

_1 


' 4 

-1 Q ■ 


-I 

B -I 

, where B = 

-1 

4 -1 


_ 0 

-I B 


_ 0 

-1 4 _ 


where h = -^^=^and the eigenvalues and eigenvectors of A are well-known. We used Householder’s 
method to tridiagonalize it. We use the same residual and orthogonality measures defined by [5] 


n = 


IIAQ-QAII2 


and 


I|/-Q^QI|2 


neP||2 


ne 
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where QAQ"^ is the computed spectral decomposition of A. 


n 

Accuracy 

measures 

Eigensolver 

QR 

CDC 

RTDC 

9 

TZ 

0.212 

0.123 

0.226 

O 

0.391 

0.334 

0.156 

25 

n 

0.202 

0.226 

0.224 

o 

0.292 

0.411 

0.174 

100 

7^ 

0.260 

0.199 

0.190 

o 

0.269 

0.163 

0.113 

400 

n 

0.235 

0.177 

0.318 

o 

0.081 

0.075 

0.068 


Table 7.1 


Accuracy measures of each algorithms 


We now compare computational complexity. Since we have to use QR (or similar) method to 
solve the subproblems, we need O(n^) operations. In this experiment, we used implicit QR for all 
algorithms. Since we can divide A into the three roughly equal sized sub-matrices, RTDC takes 
I flops of CDC to calculate the eigenvalues. For eigenvectors, CDC costs + 0{n^) arithmetic 
operations. However, RTDC takes + 0[v?) operations when we do not fix the orthogonality 
problem or use the first method in Section 16.11 If we fix the orthogonality problem by the second 
method, we need extra operations. However, since eigenvector calculations consist of nothing else 
but matrix multiplications in (12.21) . it can be effectively parallelized. Since our algorithm can calculate 
the eigenvalue twice as fast as the original CDC, our algorithm is easy to parallelize. 

8 . Conclusion. We have extended the work of Cuppen’s divide and conquer algorithm to rank 
two modification. Unlike the CDC, the deflation steps in our algorithm are nontrivial and the number 
of the eigenvalues vary in each interval li. We have shown how to classify and deflate the case when 
eigenvalues of Ti coincide with those of A. Also, we have found an algorithm to find the eigenvalue 
distribution. In the original CDC algorithm, each subinterval contains exactly one eigenvalue, but 
in our algorithm such relation no longer holds. Using the secular equation we have shown that each 
interval has either no eigenvalues, one or two eigenvalues. Eigenvectors can also be computed by the 
algorithm. However, to improve the orthogonality problem, we suggested to apply the idea used in 
the rank one modification. Since the orthogonality problem does not arise in this case, we can apply 
this algorithm in a recursive way. We believe that our algorithm can be effectively parallelized, and 
we will leave this as a future work. 
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